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ABSTRACT 

A general multi-block three-dimensional volume grid generator is presented which is suitable for Multi- 
Disciplinary Design Optimization. The code is timely, robust, highly automated, and written in ANSI “C" for 
platform independence. Algebraic techniques are used to generate and/or modify block face and volume 
grids to reflect geometric changes resulting from design optimization. Volume grids are generated/ modified 
in a batch environment and controlled via an ASCII user input deck. This allows the code to be 
incorporated directly into the design loop. Generated volume grids are presented for a High Speed Civil 
Transport (HSCT) Wing! Body geometry as well a complex HSCT configuration including horizontal and 
vertical tails, engine nacelles and pylons, and canard surfaces. 

INTRODUCTION 

The continuing increase in computer speeds and the advancement of numerical algorithms has helped to 
increase the interest in Multi-disciplinary Design Optimization (MDO). MDO is a methodology for the 
design of complex engineering systems and subsystems that coherently exploits the synergism of mutually 
interacting phenomena. The process involves analyzing independent solutions generated after perturbing the 
design variables of a baseline system. 

Typically, the MDO process proceeds in an iterative manner. Each cycle, at a minimum, includes the 
generation of numerical solutions, determination of design sensitivities, and system optimization. It is 
apparent that a large number of cycles may be necessary to complete the optimization, and as a result, the 
process must be capable of efficient execution in the batch environment with as little human intervention as 
possible. 

Perturbations are based on design sensitivities which are derivatives of design parameters (e.g. lift) with 
respect to the design variables (e.g. airfoil camber). In practice, design sensitivities from the various 
disciplines used in the system analysis are obtained separately. These are then combined using the chain rule 
to obtain the global sensitivity of the system. The derivatives may be obtained analytically, by differentiating 
the analysis code(s), or numerically, using finite differences. When using finite differences, however, it is 
often difficult to determine the appropriate step size to be used for a given system. 

In many applications of MDO, computational fluid dynamics (CFD) is an integral part of the design 
process. As a result, it is necessary to possess a rapid and highly automated grid generation capability which 
produces changes in the surface and volumetric grids to reflect the perturbations of the baseline system. It is 
also desirable to make grid modifications within the design cycle. Many design parameters, such as lift and 
drag, are fairly well established after a small number of flow iterations. With geometric changes occurring 
within the design loop, partially converged flow solutions can be used as a starting solution for the next cycle 
thereby reducing the overall computational effort 

Current grid generation techniques have been strongly shaped by a push to develop interactive tools 
which aid in the discretization of computational domains. While most of these tools are well suited for the 
generation of grids about unique configurations, their generality requires a large degree of human 
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interaction. Even for the smallest changes, like those resulting between MDO cycles, a great deal of input of 
varying degree is needed to redefine the computational structure. Recently, some tools have incorporated 
useful parametric capabilities, but continue to rely on interactive software control. These methods fall short 
of satisfying the rapid, "hands-off' grid generation needs of CFD in the MDO cycle. . 

In this paper, the development of a computer program specifically designed for MDO grid generation is 
described. The code is capable of generating CFD quality surface and volume grids, as well as analytic gnd 
sensitivities with respect to the design variables. In addition, existing grids generated with other gnd 
generation systems may be modified to reflect the geometric changes defined in an automated design and 
optimization cycle. This program has been applied in the generation and modification of both inviscid and 

In the sections to follow, the approach taken to develop the Coordinate and Sensitivity Calculator for 
Multi-disciplinary Design Optimization (CSCMDO) computer program is described Also discussed are 
results from aerospace test cases involving High Speed Civil Transport (HSCT) configurations and uture 
directions of code development. 


APPROACH 

The CSCMDO code includes Baseline 

specialized features required by MDO Surfaced) 

grid modifications. These features are 
added to compliment standard volume 
grid generation methods. The code is 
controlled via an ASCII user input file 
for execution in a batch environment. 

The code developed herein is capable of 
modifying any of the six faces of a 
block, representing a computational 
cube, to reflect geometric changes in the 
optimized system as defined by input 
surface(s). 

Figure 1 shows the layout of a 
typical MDO design loop. Information 
input from outside the loop is generated 

one time before the loop is initiated. Figure 1 Integration of grid generation into the design loop 

This information includes the baseline B 6 66 

geometry, baseline CFD grid, and a user __ 

input file. The design loop is totally self contained and therefore requires no human intervention. 1 he 
CSCMDO code operates within the loop to provide automated volume grid generation / modification within 
the each dcsi&n cycle. 

The baseline surfaces are provided in the form of a structured mesh of discrete point data. The number 
and distribution of points defining the surfaces are not required to match those of the desired CFD gnd. 
However, sufficient point resolution must be provided so as to adequately define all surface curvature. Some 
methods of surface modification described herein will require that the surface mesh topology match that ot 
the desired CFD grid. Also, points of interest, such as a wing crank location, must fall along isoparametric 
lines in the surface definition. Point continuous surface/surface intersections are desirable. In the event that 
point continuity cannot be provided at surface/surface intersections, intersection curves on each surface 
should have sufficient resolution such that interpolation errors resulting from the inconsistent data do not 

contaminate the quality of the resulting CFD grid. ... .. „ . . . , 

The baseline grid may be in the form of a baseline volume gnd or, in some instances, discrete gnd point 
data for the six faces of each block. The original grids can be generated using any structured grid 
generation package. File formats used are widely used throughout the field of CFD. . 

Within the design loop, the optimization process is to provide modified surface geometry definition, 
again as discrete point data 4 ' 6 . The changes represented by the modified surfaces are assumed to be small 
so as not to violate the original topology definition. In the event the modifications do violate the topology, 
grid quality checks provide return codes to the software controlling the loop for appropriate action. 

Standard algebraic two-dimensional(2D) and three-dimensional(3D) grid generation functionality is 
implemented for the generation of faces and blocks. Blocks may be automatically subdivided using 
intermediate face "hard planes" to increase the quality of 3D interpolation. Hard planes may be input or 
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calculated internally. 
FACF. MODIFICATIONS 


Face modifications are provided using a variety of methods to be described in the following sections. 
These methods include simple input of a modified face, parametric updates to a modified surface, projection 
to a modified surface, and deformations conforming to a modified surface. 

PARAMETRIC: 

The parametric updates 7 include methods utilizing mappings to uniform parameter space (UPS), as well 
as arclength parameter space (APS). APS is provided in the event that input surface parameterization is not 
smooth and orthogonal. The use of UPS with such a surface could produce an unacceptable CFD grid. A 
detailed explanation is contained in ref. 7. APS seeks to alleviate such shortcomings. 

These methods use a forward mapping into the selected parameter space where the grid generation takes 
place. Parameter values are then mapped back to physical space according to the modified surface. 

The mapping to and from both UPS and APS, as well as a description of the parameter space grid 
generation technique of reference 7, are briefly described here for completeness. 

Forward mapping involves the transfer of the surface grid from physical space = {X,Y,Z} T j to a 

parameter space = {(/,V'}^'j. The surface is a parametric surface such as a bilinear, bicubic, or NURBS 

representation. UPS is the primary type of parameter space used in grid generation, and is simply 
constructed using the surface grid indices according to 

Uij=i„ V,j = j. ( 1 ) 

APS is constructed using the surface arclengths in each computational direction as 

Ui.j = Vi-\.j +e i.j' 

V i,j = V i-U + f‘.J' 

where 

e i,j = |^i./ 

/ i.j ~ | ^i.j 

Grid generation is the same for both parameter spaces. The edges of the CFD face are determined in the 
appropriate parameter space. In CSCMDO, the orientation of the CFD face to the background surface is 
automatically determined. A manual override is provided in the event the orientation search fails. Once the 
edges are known, the interior of the face is generated using arclength based transfinite interpolation (TFI). 
With TFI, the interior grid points are defined by the edge and comer contributions according to 

Wl.m = + b lm W lm2 + C lm W, Um + d, m W l2m 

( 4 ) 

~ ^l,m^l.m W l2,n2- 

The subscripts (11,12) and (ml, m2) are the minimum and maximum indices of the face in the u and v 
directions respectively. The variables a,b,c and d are the blending functions subject to the following 
conditions. 


( 2 ) 


(3) 


a l,ml ~ C ll,m ~ ty'„2 ~ d/ 2>m ~ 1. 
a l.m2 = c l2,m = 1 = ^il.m = 0. 


(5) 
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A very robust set of blending functions has been proposed by Soni^. These blending functions are based on 
the normalized arclength of the edges in physical space, and are defined as 


where 


a l,m ~ Vl,m * 
C l,m ~ 1 “ 


bl,m “ ♦ 

( 6 ) 

^l,m £/,m» 


e £/,ml £f,ml } 




hn.m + 8t,mi \hn.m ~ ^n,«} 

m.rn = 1 L 


p 






( 7 ) 


The values of g[ m and /i/ >m are the normalized edge arclengths in physical space. With this interpolation 
complete the grid can be mapped back from the parameter space to the physical space. 


Backward mapping is used to transfer the grid points from the parametric space = {u,. m ,v /m } T j to 

the physical space = {.*, m ,y im ,z /(B } r ^. The backward mapping involves three steps. The first step is a 
search of the surface parameter space to determine the local surface patch which includes the current grid 
point. The second step is to compute the local patch coordinates (si n ,t lm ) for the coordinates in parametric 

space (w, *). For UPS this is relatively simple, where 


h,m ~ ^ l,m ~ ^ i,j • 


( 8 ) 


For APS, the local coordinates are computed using a set of quasilinear equations. The UPS parameter values 
of U t i and V t j are known for each comer of the patch and the patch grid point is known in the 

interior of the patch. With this information a system of two equations with two unknowns is formed, where 


U l,n + S l.m (l ~ 

+ (l — s l,m ) t l,mUi,j + 1 + s l,mh,m^i+\,j+\ > 

( 9 ) 

v /,m = (l ~ S l,m )(l - h,m )^i.j + - 

■*■(1 — s l,m • 

The solution to these equations is obtained using the Newton-Raphson method and gives the local 
coordinates of the patch The final step is to compute the physical coordinates (r, m ) of the grid 

point on the local surface patch using bilinear or bicubic interpolation. 

PROJECTION: 

Projection of a face allows for the generation of a CFD grid over several surfaces, while also avoiding 
problems associated with surface parameterization mentioned in the previous section. The projection 
algorithm^ used in CSCMDO is described briefly for completeness. Again the surface is approximated by a 
parametric surface as, 
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( 10 ) 


R(w)={x(w),y(w),z(w)} T , 
vv = {u.v} 7 . 

The process of projecting a point, r, onto a surface /?(w)> involves finding w such that the distance, d, 
between f and R(w) is minimized such that it does not violate the limits of w. The distance, d, can be 
written in terms of parameter w as, 


d 2 (w) = f(w) = \R(w)-r\\ (11) 

To minimize d, Eq. (11) must be minimized with respect to w. This is accomplished by setting the gradient 
of /. V/(w?), equal to zero, as 

V/(w) = G,(w) = -^ = 0, 

' ( 12 ) 

Solution of the above nonlinear system of equations for w results in the projected point. 

As with the parametric update, a surface is represented as a set of bilinear or bicubic patches. For each 

patch, R(u,v), is approximated in terms of its parameters as, 


R («. v) = (1 - «)(1 ~ v)Rij + u(l - v)R mj 
+(1 -u)vR ij+l +uvR MJ+l . 


(13) 


Combining Eqs. (12) and (13) yields the following system of nonlinear equations, 


where. 


-F} = 0, ^-{R(w)-r} = 0, (14) 


^ = (1 ■ - v)(* i+u - Rij ) + v(/? 1+1>>+1 - R iJ+l ), 

^ = (1 -upi, /+1 -Rij) + u(R, +u+i -R MtJ ). 


(15) 


Again the system is solved by the Newton-Raphson method. An additional method of projection available in 
CSCMDO is the projection onto an original surface(s) and evaluation based on a new surface(s). This 
method obtains the parametric coordinates by projecting each face grid point onto the baseline surface. The 
parametric coordinates are then mapped to the physical space using the surface(s) driving the modification. 
The resulting grid maintains the same characteristics of the original grid. The shortcoming of this method, 
however, is die required existence of both the original and modified surfaces at the time of projection. 

DEFORMATION: 


Deformation of a face is key in providing reusability of an existing CFD grid. When a face is modified, 
a comparison to the original face edge points is made according to. 


NR - R mo & ifitd ~ ^original > ( 16 ) 
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Rearranging Eq. (16), the modified edge is given by. 


^moiifitd ^ original + AR . 


(17) 



a) Original smoothed, orthogonal face 


If all edges have been modified, Eqs. (16) and (17) are applied to the interior grid points to yield the 

modified face. In the event that less than 4 edges of a face are modified, unmodified edge AR values in Eq. 

(17) are obtained by blending the ending values of 

A R along the edge. The edge AR values are then 
interpolated into the interior for use in Eq. (17). This 
method maintains the grid quality characteristics of 
the original CFD grid and provides time saving 
reusability of an optimal CFD grid. 

Figure 2 shows the method applied to a C-H grid 
around a circle enclosed in a rectangular outer 
domain. The original grid in figure 2-a shows 
smoothness and orthogonality characteristics obtained 
through partial differential equation (PDE) 
smoothing. Figures 2-b and 2-c show a modification 
to the circular inner edge. In Figure 2-b the edge is 
modified and the face is re-initialized with arclength 
based TFI. Figure 2-c shows the effect of a face 
deformation governed by the same change to the 
circular edge. Note the preservation of the 
smoothness and orthogonality characteristics in figure 
2-c. This method is extended to 3D as will be 
discussed in later sections. 

BLOCK MODIFICATIONS 

Volume modifications are accomplished using 
standard algebraic methods for the re-initialization of 
the block interior, or by deformation of the original 
block interior based on changes defined by the six 
faces. 

GRID GENERATION: 

The generation of the block interior is conducted 
in a manner which is an extension of the 2D methods 
of arclength based TFI described above. As 
previously mentioned, a block may be automatically 
subdivided by inserting intermediate faces, or “hard 
planes” within the block. Each hard plane acts to 
subdivide the block providing added control over the 
quality of developing interior. Hard planes may be 
directly input, as in the case where a PDE smoother 
was used to obtain an optimal grid, or generated 
internally using the 2D methods described above. 

DEFORMATION: 

A method similar to that described for face 
deformation is extended to 3D for use in the 
deformation of entire blocks. Once the deformations 
are calculated for each of the six faces of a block, 
arclength based TFI is used to blend the deformations 



b) Modified face with algebraic TFI 



c) Modified face with face deformations 
Figure 2 Example of Face deformation 
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into the block interior. Deformations are then applied to the original grid points to obtain a deformed block. 
Again this method of deformation facilitates the reusability of an existing volume grid. The overall quality 
and characteristics of an existing block are maintained after face modification in the modified block. 

Figure 3 shows the 3D deformation applied to the same C-H topology used previously in 2D. Figure 3-a 
shows smoothness and orthogonality characteristics obtained through partial differential equation (PDE) 
smoothing for the original block. Figures 3-b and 3-c show the same modification to the cylindrical inner 
face. In Figure 3-b the face is modified with the adjacent faces and interior re-initialized with arclength 
based TFI. Figure 3-c shows the effect of face and block deformation governed by change to the cylindrical 
inner face. Again note the preservation of the smoothness and orthogonality characteristics in figure 3-c on 
the block interior. 





a) Original smooth orthogonal block b) Modified with algebraic TFI c) Modified with deformation 


Figure 3 Example of block deformation 


GRID QUALITY 

A check of grid quality measures is included in the code for diagnostic output as well as limited control when 
running in a batch environment. Cell volume and cell skewness are calculated on a block by block basis. A 
histogram of cell skewness is provided with the code output showing the percentage of the grid possessing a 
skewness over a range of 0° to 90°. Minimum cell volume, maximum cell skewness, and maximum average 
cell skewness controls are provided to trigger return codes which are used by the process controlling batch 
execution. The return codes allow for the controlling process to exit the design loop without wasting 
continued resources in the event the resulting volume grid is not suitable for use. 

SENSITIVITY 

Calculation of grid sensitivity to the design variables involves the automatic differentiation of the CSCMDO 
code. This is accomplished using the Automatic Differentiation of C source (ADIC) processor from 
Argonne National Labs. ADIC is a tool for the automatic differentiation of ANSI C programs providing 
similar capabilities for C as ADIFOR^ does for FORTRAN77. The ADIC tool operates on a given source 
code using a specified set of dependent and independent variables, to produce an augmented C source that 
computes not only the original scalar result, but also the partial derivatives of all of the specified dependent 
variables with respect to the independent variables. ADIC employs the Sage++ programming environment 
developed at Indiana University. 

The cost of generating the sensitivities with the ADIC processed version of CSCMDO (CSCMDO. AD) has 
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proven to be negligible when compared to the overall cost of a CFD design cycle. With an ADIC processed 
source code, the ratio between the cost of evaluating a gradient with n components and the cost of evaluating 
the original underlying scalar function is not n+1, but is bounded by 5 at most^. Preliminary testing of 
CSCMDO.AD has produced reasonable results. These results have been compared to finite difference 
calculations and show favorable agreement. 

VISCOUS GRID CAPABILITIES 

The methods presented herein have been successfully applied to viscous grids. Care must be taken to provide 
surface geometry containing sufficient resolution to capture surface curvature. The use of block 
deformation is very useful in the reuse of viscous grids which were originally refined using PDE techniques. 

RESULTS 

CSCMDO has been extensively tested on aerospace configurations at the NASA Langley Research 
Center’s (LaRC) Geometry Laboratory (GEOLAB). The test cases range from simple wing/body 
configurations to full HSCT geometry with tail surfaces, engine nacelles, and canards. The code is also used 
outside of the design loop in the GEOLAB for the rapid modification and quality check of existing CFD 
volume grids. 

Execution time is difficult to determine for the code due to the wide variety of execution schedules 
available to the user. Typical timing, however, can be given as an example of two types of configurations. 

The first is that of a simple double delta wing/body aircraft. The computational domain was broken into 
two C-H block with a total of 350,000 grid points. Seven modified geometries were supplied along with the 
baseline and the block topology. The geometry modifications included changes to wing sweep, camber, 
span, crank location, etc.. All eight volume grids, baseline and seven modifications, were generated 
sequentially on a Silicon Graphics Onyx workstation in double precision in under three minutes. 



Figure 4 Complex HSCT configuration 

The second case is that shown in figure 4. The configuration includes a double delta wing, tail surfaces, 
two flow through engine nacelles with a plug insert, engine pylons, and a canard. The computational 
domain contained seventeen blocks with at total of over 1.5 million grid points. This complex case required 
on the order of three minutes for the generation of a single complete volume grid with modifications to a 
variety of surfaces. 

Figure 5 presents the results of the case involving the complex HSCT in figure 4. Figure 5-a shows every 
third grid line in each computational direction on the surface of the configuration. A cut from the 3D 
volume is shown in figure 5-b. The view is looking downstream at the nose of the aircraft. 

FUTURE DIRECTIONS 

Future plans include the addition of PDE smoothing of block faces and entire blocks within the CSCMDO 
code. Also to be included are capabilities for complete domain decomposition within the code. Presently the 
block topology, number of points, and relative spacing are required input to CSCMDO. Future plans are to 
incorporate methods of defining these key features of a computational domain within the program. 
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CONCLUSIONS 


A CFD volume grid generator has been presented for which is designed around the needs of multi- 
disciplinary design optimization. The code employs methods which make it ideal for the type of grid 
modifications required by MDO. The code is robust, timely, platform independent, and ready to meet the 
needs of multi-disciplinary design. 
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b) Grid from cut normal to flow 
Figure 5 Examples of complex volume grid 
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